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1.  Objectives 

The  main  focus  of  the  work  was  the  development  of  a  fast  time-domain  algorithm  for  solv¬ 
ing  Maxwell’s  equations  in  the  integral  form  for  general  dispersive  media  with  both  penetrable 
and  impenetrable  surface  and  bulk  volumetric  properties.  The  main  motivation  for  this  work 
was  the  need  for  an  integral  equation  solver  applicable  to  the  description  of  propagation 
of  wide-band  pulses  through  different  types  of  dispersive  media  including  those  which  are 
of  interest  in  the  context  of  foliage  penetration. 

During  the  course  of  the  project  we  have  been  working  closely  on  related  applications 
of  the  fast  integral  equation  solver  with  the  Brooks  Air  Force  Base  Mathematical  Products 
Group  lead  by  Dr.  R.  Albanese. 


2.  Main  accomplishments 

We  summarize  the  main  developments  carried  out  on  this  project. 

We  constructed  a  new  time-domain  algorithm  which  is  particularly  well  suited  to  gen¬ 
eral  dispersive  media  with  both  penetrable  and  impenetrable  surface  and  bulk  volumetric 
properties.  The  algorithm  is  described  in  [5]  and  in  the  first  yearly  report.  The  main 
features  of  the  algorithm  are  as  follows: 

(i)  The  algorithm  complexity  does  not  depend  on  the  degree  of  dispersion.  Its  small  com¬ 
putational  cost  ( 0(Nt  Ns  log  Nt  log Ns)  and  0(Nt  log  Nt  log  Ns)  for  volume  and 
surface  problems  respectively,  where  Nt  and  Ns  denote  the  number  of  temporal  and 
spatial  samples)  is  achieved  through  the  simultaneous  application  of  Fast  Fourier  Trans¬ 
forms  in  space  and  time.  We  stress  that  algorithm  is  applicable  to  both 

(il)  conventional  volumetric  integral  (involving  free  space  Green’s  function  and  volu¬ 
metric  discretization)  equations  for  dispersive  media 

and  to 

02)  (  much  more  difficult  in  terms  of  formulation,  analytical  and  numerical  imple¬ 
mentation,  but  significantly  more  beneficial  in  terms  of  algorithm  complexity) 
equivalent  current  surface  integral  equations  for  a  general  piecewise  homogeneous 
dispersive  dielectric  medium  with  dispersive  thin  sheet  material  interfaces.  (We 
note  that  the  main  part  of  our  effort  was  devoted  to  such  integral  equations.) 

(ii)  The  algorithm  is  based  on  a  new  formulation  of  integral  equations  specially  tailored  to 
problems  involving  dispersive  media.  The  formulation  is  both  general  and  significantly 
simpler  than  the  conventional  approaches:  instead  of  using  the  customary  integral 
equation  operators  involving  the  Green  function  and  its  derivatives,  we  construct 
supplemental  integral  equation  operators  equal  (a)  to  the  Fourier  transform  of 
the  dispersive  medium  Green  function,  (b)  to  the  Fourier  transform  of  the  product  of 
the  dispersive  medium  Green  function  with  the  frequency  dependent  dielectric  per¬ 
mittivity,  and,  (c)  to  the  Fourier  transform  of  the  product  of  the  dispersive  medium 
Green  function  with  the  inverse  of  the  dielectric  permittivity.  An  important  bene¬ 
fit  of  such  an  approach  is  that  the  resulting  integrals  involve  only  single  (and  not 
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double)  time  convolutions.  The  formulation  is  applicable  to  systems  involving  bulk 
dispersive  regions  and  thin  dispersive  sheets  represented  as  interfaces.  We  have  car¬ 
ried  out  complete  analytical  calculations  and  corresponding  numerical  procedures  for 
the  evaluation  of  matrix  elements  of  the  integral  operators,  executed  in  the  frame¬ 
work  of  the  full  Galerkin  scheme  in  space  and  time  variables,  for  the  “conductive 
Debye  medium”  (i.e.,  for  a  medium  with  the  electric  permittivity  given  by  the  Debye 
formula  supplemented  with  a  term  responsible  for  the  medium  conductivity).  The 
procedure  employs  a  suitable  contour  integration  around  singularities  of  the  effective 
Green  function  operators  in  the  complex  frequency  plane. 

(Hi)  The  algorithm  does  not  require  analytical  representation  of  the  (dispersive)  medium  Green 
function  and,  in  particular,  can  be  used  with  a  medium  Green  function  given  in  a  tabular 
form. 

(iv)  Our  formulation  provides  a  framework  and  a  practical  tool  for  an  accurate  numerical 
simulation  framework  for  a  variety  of  problems  involving  wide  band  pulse  propaga¬ 
tion,  including  wide  band  antenna  pattern  simulation  or  propagation  of  narrow  pulses 
through  dispersive  media.  Such  applications  require  the  numerical  capability  of  prop¬ 
agating  an  ultrawideband  pulse  in  a  dispersive  environment  over  several  absorption 
lengths.  In  particular,  such  numerical  capability  appears  to  be  necessary  in  studies 
of  an  intriguing  phenomenon  of  a  Brillouin  precursor,  a  propagating  structure  which, 
because  of  its  non-exponential  peak  decay,  may  lead  to  multiple  potential  applications, 
including  detection  of  objects  in  cluttered  environments. 

We  note  that  our  work  on  this  project  had  a  significant  impact  on  the  CEM  academic 
and  industrial  community.  In  particular  an  early  version  of  the  block  Toeplitz  time-domain 
integral  equation  solver  we  proposed  has  been  independently  verified  and  adopted  as  the 
Time  Domain  Adaptive  Integral  Equation  Method  by  the  University  of  Illinois  researches 
who  constructed  a  code  based  on  this  formulation  [3].  Also  researches  in  Motorola  and 
Cadence  [4]  apply  the  method  in  the  context  of  electronic  packaging. 
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3.  Technical  details 


In  this  Section  we  describe  in  more  detail  the  relevant  aspects  associated  with  the 
formulation  and  implementation  of  the  approach. 


(A.)  Time-domain  electric-field  integral  equations  for  dispersive  media 

The  integral  equation  formulation  we  developed  is  particularly  suitable  to  problems  involving 
dispersive  media.  We  constructed  a  general  set  of  integral  equations,  significantly  simpler  than 


the  conventional  formulation ;  instead  of  using  the  customary  integral  equation  operators 
equal  to  the  Green  functions  g(oj,  r)  and  their  derivatives,  we  use  supplemental  integral 
equation  operators  K0(t,  r),  Kft,  r),  and  K2(t,  r)  equal  to  the  Fourier  transforms  of  the 

Green  function  of  the  dispersive  medium  and  their  products 
dielectric  permittivity 

with  frequency  dependent 

g{t,r)  =  J  ^~e_1Wt5(w,r)  , 

(3.1a) 

9'(i'r)  =  /  2,Z(u,)e 

(3.1b) 

52  M  =  J  ^I(uj)e-luJtg(uj,r)  . 

(3.1c) 

An  important  benefit  of  such  formulation  is  its  generality,  simplicity  and  efficiency.  The 
resulting  integrals  involve  only  single  (not  double)  time  convolutions,  i.e.,  the  resulting 
integral  equations  have  the  schematic  structure  J  K(t  —  t')X(t')dt'  —  B(t),  instead  of 
fK1{t-t")K2(t"-t,)X(t')dt"dt'  =  B(t)  . 

In  terms  of  g(t,  r),  gft,  r),  and  g2(t,  r)  the  surface  integral  equations  for  a  general 
piecewise  homogeneous  dispersive  dielectric  medium  with  dispersive  thin  sheet  material 
interfaces  take  the  form 


Udt' 

+  Jdt'  J d2S' 


R(t  —  t',r)  U(t  —  t',  r)n(r)x 


U(t  —  t',  r)n(r)x  S(t  —  t’,r) 


M(t',r) 


Ku(t-t',r-r')  K12(t-t',  r  -  r') 
K21{t  -t',r  -r')  K22(t-t',  r-r') 


"  3(t',  r') 

*  E'(t,  r) " 

M(t',r') 

_Hi(t,r)_ 

(3-2) 


with  the  kernels  K-  given  by 

t 

Ku(t-t',  r-r')  =  ^  dtg(t-t' ,r  -  r')  -  V  V  •  J  drg^r  -  t',r  -  r')  , 

—  OO 


Kl2(t  —  t',  r  -  r')  =  c  1  V  x  g(t  -  t' ,  r-r')  , 

K21(t  -t',r-  r')  =  -c_1  V  x  g(t  -  t',r  -  r')  , 

t 

K22(t-t',r-T,)  =  \dtg2(t-tl,r-r,)-ixX-  f  dr  g(r  -  t',r  -  r’)  ,  (3.3) 

c  p  j 

—  OO 


3 


In  the  above  equation  R(t,  r),  S(t,  r),  and  U(t,r )  represent  the  Fourier  transforms  of  the 
respective  functions  R(u,  r),  S(u,  r),  and  U(oj.  r)  representing  the  dimensionless  electric, 
magnetic,  and  “cross-”  resistivities  of  the  mathematically  infinitely  thin  sheet  material 
interfaces  [1]  between  two  homogeneous  volumetric  regions,  and  J(t,r)  and  M(t,  r)  are  the 
unknown  vector  functions:  equivalent  electric  and  magnetic  currents  on  material  interfaces. 

We  note  that  in  (3.3)  we  assumed,  for  simplicity,  that  =  const.  Generalization 
to  dispersive  ]l(u)  is  straightforward. 

An  alternative  set  of  equations,  obtained  by  taking  the  negative  time  derivative  of 
(3.2)  is 


dtR(t-t',r) 
dtU(t  -  t',  r)  n(r)x 


dtU(t  -  t',  r)  n(r)  x 
dtS{t-t\  r) 


Kn(t  -  t1  ,r  -  t') 


K2i(t-  t’,r-  r') 


K12(t-t',r-r') 

K22(t-t',r-r') 


J(t',r) 

M(i',r)_ 


Ei(t,  r) 

H^r) 


,  (3-4) 


with 


Kn(t  —  t',r  —  r')  =  ^  d\  g(t  —  t',  r  —  r')  —  V  V  •  gx(t  -  t',r  -  r')  , 

Ku(t  -  i',  r  -  r')  =  c_1  V  x  dt  g(t  r')  , 

K2l{t  -t',r-  r')  =  -c_1  Vx9(  g(t  -t\  r  -  r')  , 

K22(t  -t',r-  r')  =  -^d?g2(t  -t’,r  -  r')  -=VV-j(t-  t\ r  -  r')  .  (3.5) 

c  i  L 

We  stress  that  both  the  fast  FFT  solution  scheme  [5]  and  the  complete  Galerkin 
discretization  method  in  space  and  time  which  we  discuss  below,  are  applicable  to  Eqs.  (3.2) 
as  well  as  (3.4). 
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(B.)  Complete  Galerkin  discretization  procedure  in  spatial  and  temporal  vari¬ 
ables 

We  developed  an  accurate  method  of  the  treatment  of  the  causal  behavior  of  the  time-domain 
Green  function  achieved  through  the  utilization  of  the  complete  Galerkin  procedure  in  spatial 
and  temporal  variables  in  the  evaluation  of  impedance  matrix  elements. 

It  has  been  known,  in  general,  that  the  use  of  the  Galerkin  method  constitutes  a  key 
element  in  achieving  numerical  stability  of  the  solution  schemes  for  time-domain  integral 
equation  approaches.  Although  for  one-dimensional  wave  propagation  problems  implemen¬ 
tation  of  the  causality  of  the  Green  function  in  the  context  of  the  Galerkin  discretization 
is  relatively  straightforward,  implementation  of  the  causality  in  three  dimensions  presents 
a  much  more  complicated  task  for  the  following  two  reasons: 

(a.)  The  Green  function  depends  on  the  relative  time  r  and  distance  r  —  |ri  —  r2 1-  If 
the  basis  functions  have  a  certain  spatial  extent,  the  signal  from  one  basis  function 
may,  in  general,  arrive  only  to  a  portion  of  the  other  basis  function.  Implementation 
of  causality  requires  careful  implementation  of  all  relevant  radial  integrals  over  the 
spatial  coordinates  iq  and  r2  associated  with  the  basis  functions  entering  the  matrix 
element. 

(b.)  The  three-dimensional  Green  function  has  a  singular  term  associated  with  1/r  = 
1/ |r i  —  r 2 1 ,  which  is  in  practice  rather  difficult  to  take  into  account  in  the  context  of 
three-dimensional  integrals  over  vector  variables  iq  and  r2. 

We  implemented  routines  for  an  efficient  and  accurate  evaluation  of  the  spatial  separation 
densities  ( Pap{r ))  (defined  as  angular  averages  of  convolutions  of  pairs  of  spatial  basis 
functions)  for  Rao-Glisson-Wilton  basis  functions.  In  Figs,  la  through  If  below  we  display 
some  numerical  results  for  the  separation  densities  calculated  according  to  our  prescription, 
associated  with  various  components  of  the  Green  function.  Numerical  values  of  the  sepa¬ 
ration  densities  obtained  using  our  approach  compare  almost  exactly  with  the  analytically 
equivalent  but  numerically  much  more  intensive  expressions  involving  direct  3-dimensional 
numerical  integration  quadratures. 
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Fig.  la:  Numerical  values  of  separation 
densities  associated  with  9  components 
of  the  tensor  Green  function  for  the  con¬ 
figuration  of  2  triangles  in  the  far-held 
zone. 
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Fig.  lb:  Numerical  values  of  separation 
densities  associated  with  9  components 
of  the  vector  Green  function  for  the  con¬ 
figuration  of  2  triangles  in  the  far-held 
zone. 
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Fig.  lc:  Numerical  values  of  separation 
densities  associated  with  9  components 
of  the  tensor  Green  function  for  configu¬ 
ration  of  2  adjacent  triangles  (near  held 
zone). 
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Fig.  Id:  Numerical  values  of  separa¬ 
tion  densities  associated  with  9  compo¬ 
nents  of  the  vector  Green  function  for 
the  conhguration  of  2  adjacent  triangles 
(near-held  zone). 
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Fig.  le:  Numerical  values  of  separa¬ 
tion  densities  associated  with  9  compo¬ 
nents  of  the  tensor  Green  function  for 
the  configuration  of  2  overlapping  trian¬ 
gles  (self-field  zone). 
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Fig.  If:  Numerical  values  of  separa¬ 
tion  densities  associated  with  9  compo¬ 
nents  of  the  vector  Green  function  for 
the  configuration  of  2  overlapping  trian¬ 
gles  (self-field  zone). 


More  details  on  this  part  of  the  project  were  reported  in  the  second  yearly  progress  report 
and  published  in  [5]. 
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(C.)  Formulation  for  a  conductive  Debye  medium 

We  have  carried  out  complete  analytical  calculations  for  a  specific,  and  general  model  of 
a  “conductive  Debye  medium” ,  with  the  electric  permittivity  given  by  the  Debye  formula 
supplemented  with  a  term  responsible  for  conductivity, 


«(<*>)  =  Coo  + 


1  —  ITU 


.  V 

+  1-  . 
u 


(3.6) 


In  Eq.(3.6)  es  and  are  the  static  and  optical  relative  permittivities,  r  is  the  relaxation 
time,  and  v  =  a/ e0  =  cZ0  a  is  a  parameter  proportional  to  the  conductivity  a,  having  the 
dimension  of  frequency.  Analyticity  requires  real  es  and  e^,  and  we  assume,  additionally, 
that  es  >  e^. 

The  parameterization  (3.6)  reduces  to  the  following  cases: 


(non-dispersive  medium)  for  u  —  0,  r  =  oo  , 


e(u)  =  { 


6„  +  — — r-22-  (Debye  medium)  for  v  =  0  , 

1  —  ITU 
V 


(3.7) 


u 


(lossy  medium)  for  r  =  oo  . 


To  exhibit  the  analytic  structure  of  the  electric  permittivity  function  we  now  rewrite 
Eq.(3.6)  as 

,  (w  -  wx)  (u  -  u2) 

-  €oo 


U  (u  —  u3 ) 


where 


(3.8) 


u n  = 


2  eoo  r 


UJn  — 


2eoo'r 


Ur,  = - . 

3  T 


~(es  +  vr)  +  V(es  -  VT)2  +  4  (es  -  eTO)  vr 
-(es  +  vt )  -  \/(es  -  VT)2  +  4  (cs  -  O  vt 


(3.9a) 

(3.9b) 

(3.9c) 


Thus,  as  required  by  the  analyticity  (causality)  properties,  e(u)  is  analytic  in  the  upper 
half-plane  of  u;  it  has  two  poles  (at  zero  and  on  the  negative  imaginary  axis,  at  u3),  as 
well  as  two  zeros,  also  located  on  the  negative  imaginary  axis. 

The  wave  number  is  then 


k(u) 


u 


\J(.{u)  IJ,(u) 

C  V 


1  u  (u  —  uq)  (u  —  u2) 


U  —  Uo 


(3.10) 


where  c  =  1/ ^/eo  Ho  and 


(3.11) 
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The  expression  (3.10)  has  thus  branch  points  at  u  =  0  and  u  =  wn,  n  =  1,  2, 3. 

Since,  for  typical  parameters,  we  have 

0  <  l^xl  —  1^3 1  —  1^2 1  >  (3-12) 

it  is  convenient  to  define  branch  cuts  in  the  u  plane  as  in  Fig.  2. 

We  have  constructed  both  analytical  representations  and  numerical  procedures  for 
computing  matrix  elements  of  the  integral  operators  for  the  above  mentioned  medium, 
employing  full  Galerkin  procedure  in  space  and  time. 

The  analytical  calculations  are  reported  in  Appendix  A1  of  the  second  yearly  report. 
The  final  result  is  a  set  of  expressions  for  matrix  elements  of  the  Green  functions  g,  gv 
and  g2  of  Eq.(3.1)  which  can  be  represented,  schematically,  as 

Kp  =  J  dti  d2G  d2i  d2f2  Tll{h)  ^a(ri)  ’  9(h  ~  *2>  ri  ~  f2)  T°(t2)  Vp(r2) 

OO 

=  J  driP«p(r))  Fir)  ■  (3.13) 

0 

Here  is  the  matrix  element  computed  between  temporal  basis  functions  T  of  time 
separation  p  At  (where  At  is  the  time  step),  and  spatial  basis  functions  >F  associated  with 
current  elements  (edges  of  the  discretization)  a  and  (3.  The  actual  matrix  elements  may 
also  involve  Green  functions  gx  and  g2,  or  their  time  derivatives  or  integrals,  and  spatial 
differential  operators  (divergence  and  curl)  acting  either  on  the  Green  functions  or  on  the 
spatial  basis  functions. 

Eq.(3.13)  represents  the  matrix  element  as  a  space  integral  of  the  “spatial  separation 
density”  (pap(r))  (expressible  in  terms  of  the  basis  functions  \FQ  and  VF^)  and  functions  /M 
constructed  from  the  Green  function  and  the  spatial  basis  functions.  All  the  dependence 
on  the  medium  properties  and  its  dispersion  is  thus  contained  in  the  functions  /ffi 

In  Appendix  A1  of  the  second  yearly  report  we  described  in  detail  evaluation  of  the 
spatial  separation  densities  ( p )  and  provide  semi-analytic  expression  for  their  efficient  com¬ 
putation.  We  also  give  a  details  of  the  procedure  for  computing  the  functions  /ffi  in  the 
general  case  of  the  conductive  Debye  medium  (Eq.(3.6))  they  are  given  as  frequency  in¬ 
tegrals  over  the  frequency-domain  Green  functions  and  Fourier  transforms  of  “temporal 
separation  densities”  7 ?M(t),  defined  in  analogy  to  (pag)  as  convolutions  of  temporal  basis 
functions.  The  frequency  integrals  are  evaluated  by  deforming  the  integration  contour  to 
one  encircling  singularities  of  ?(w)  and  which  are  located  on  the  imaginary  u  axis; 

hence  the  expressions  involve  also  Laplace  transforms  of  the  temporal  separation  densi¬ 
ties  rj'l(t),  which  can  be  evaluated  analytically.  The  resulting  frequency  integrals  can  be 
efficiently  and  accurately  computed  as  numerical  quadratures.  Some  of  these  integrals, 
however,  have  to  be  treated  with  special  care:  this  applies  to  cases  when  a  pole  in  the 
integrand  coincides  with  a  branch  point  (which  is  also  an  essential  singularity). 
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Fig.  2:  The  singularity  structure  of  the  electric  permittivity  function  and  the  integration 
used  in  evaluation  of  the  Green  function. 


We  present  here  more  details  of  the  computation  of  the  time-domain  Green  functions 
g(t,  r),  gi(t,r)  and  52(t,  r)  of  Eqs.  (3.1)  in  the  case  of  the  conductive  Debye  medium 
(Eq.(3.6)).  First  we  consider  the  function  g(t,  r)  of  Eq.(3.1a)  defined  as 

9((-r)  =  /  =  db  .  (3.14) 

It  is  convenient  to  start  with  evaluating  the  corresponding  Green  function  in  one 
spatial  dimension,  where,  in  frequency  domain 


and  thus  the  corresponding  Green  function  in  three  dimensions  is 

5(w,r)  -  -Lelfc(u,)r  =  ^-g(u,r)  . 

47rr  2-Kr  or 

Now,  for  the  one-dimensional  case,  the  integral 

do;  i 


(3.15) 


(3.16) 


g{t,r)  =  J 
iv  r 

=  2"  J 


2-k  2  k(u) 

dcu  r 


e-iwt  gifc(w)r 


W-W, 


27T  V  u  (u  —  uq)  (ui  —  u>2) 


e  lbJt  exp  i 


lu  (w  -  uq)  (w  -  Wa)  r  \  ^  ^ 


U  —  LJ  q 


is  convergent,  and  can  be  evaluated  by  deforming  the  original  integration  contour  along 
the  real  axis  on  the  w-plane  into  the  contours  C'  and  C"  of  Fig.  2.  In  the  latter  con¬ 
tours  we  include  circles  around  singular  points;  in  some  cases  these  circles  give  nonzero 
contributions. 
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By  analyzing  the  phase  changes  on  encircling  the  branch  points,  we  find  that  the 
square  root  appearing  in  the  exponential  of  Eq.(3.17)  is  real  positive  just  to  the  right  of 
the  branch  cuts,  and  real  negative  just  to  the  left  of  these  cuts.  Then,  by  changing  the 
integration  variable  to  u  —  iw,  and  denoting  uju  —  — iun  with  Reun  >  0,  n  =  1,2,3,  we 
obtain 


s  “i 

A/  .  v  _  (  r\  I  f  du  1 


u) 


e~ut  ei/(“)r/ 


u2 

'  +  / 


du  1 


27T  f(u ) 


e-ut  &if(u)r/v 


“l 

/ 


^Le-ute-if(u)r/v  _  f  JzL0-vt*-if(u)r/v 


2 7T  /(u) 


v  _  f  du 

J  2tt 


/(«) 


e  “"e 


=  v@(t~rv)  {It W) e_ut cos (/(u) v) +  t 7R e~" cos (/(u) D  •  (3-18) 


where  ©  is  the  Heaviside  step  function  and 


/(«) 


u  (u  —  u:)  (u  —  u2) 


U  —  Uo 


(3.19) 


is  defined  to  be  positive  in  the  u  intervals  occurring  in  Eq.(3.18). 

According  to  Eq.(3.16),  the  corresponding  Green  function  in  three  dimensions  is  then 


(3.20) 


Although  the  integrands  in  Eq.(3.20)  are  singular,  these  singularities  are  only  of  one-over- 
square-root  type,  and  are  thus  integrable. 

As  we  can  see  from  Eqs.  (3.3)  and  (3.5),  we  also  need  the  spatial  derivative  of  the 
function  g(t,r).  It  is  given  by 


0  Uo 


du  1 


27T  /(u) 


e  ut  cos 


(/(“>  D 
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-2^{K)  f/+/ls  755  •-<”(/<•>;) 


{/+/}se"”‘sin(/(“)D 


V 


e  ut  sin  (/(u)  ^ 


+  7— - —  ©  ft  —  — 
2nvr  \  v 


f{u)e  ut  cos  (/(u)  -)  ,  (3.21) 


where  the  integrals  are  again  convergent. 

We  consider  next  the  Green  function  ^1(t,  r)  of  Eq.(3.1b)  and  write  it  in  the  form 

»i(t,r)=  f  =  iS,/-  .  (3-22) 

J  2n  e(co)  4irr  J  2n  ue(uj) 


where  the  last  integral  is  convergent  for  u  -»  oo.  More  explicitly, 


9'(t'r)=ia‘4 ibl% 


u  —  u3 

(w  -  Wj)  (w  -  w2) 


e-iwt  gi  fc(w)r 


(3.23) 


In  this  integral  we  can  deform  the  integration  contour  as  before  and  obtain 


9i(t,r)  =  id, 


f  f  du  1  r—ut  pi  f(u)  r  /v 

2tx  p(u ) 


dyre^r  |  J  2i r  p(u) 

v  0 


—  f  _L_  P -ut  -if  (u)r/v  _  f  du  1  -ut  c-if(u)r/v 

J  27 rp(ti)e  J  27rp(u)e 


-  -dt 


Q(t  —  r/v ) 
27reoor 


27t  p(u) 


(3.24) 


where 


p(„)  =  -i  /’(„)  = 

n  u  —  u3 


(3.25) 
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The  integrand  of  Eq.(3.24)  is  singular  at  u  =  tq  and  u  —  u2,  but,  because  of  the  behavior 
of  f{u),  the  singularities  are  integrable. 

Finally,  we  consider  the  Green  function  #2(t,r)  of  Eq.  (3.1c): 

02^, r)  =  J  ^  e~'ut  eM  0(w> r)  •  (3-26) 


Since  e(co)  — >  =  const  for  u  — >  oo,  it  is  necessary  to  subtract  the  constant  term,  which 

yields  a  delta-type  contribution  to  the  Green  function.  We  have  thus 


92(t,r)=  eOQg(t,r)  +  ^  |^[e(u)-eje 


eoo  9(t,r)  + 


4nr 


r  dw 

J  27T 


(uj  -  oq)  (u  -  uj2) 


UJ  (uj  —  CJ3) 


- 1 


-iut  i/(iu>)  r/v 


,  (3.27) 


where  Eq.(3.8)  was  used  to  express  e(u>)  and  the  function  /( iuj)  is  given  by  Eq.(3.19). 
Evaluation  of  the  last  integral  presents  a  certain  difficulty,  because,  at  ui  =  u>3,  there  is 
a  coincidence  of  the  pole  in  the  integrand  with  the  essential  singularity  in  the  exponen¬ 
tial,  resulting  from  the  one-over-square-root  singularity  in  /( iw).  As  before,  based  on  the 
asymptotic  behavior  of  the  integrand  for  uj  — >  oo,  we  can  deform  in  the  integration  contour 
in  Eq.(3.27)  to  the  contours  of  Fig.  2.  In  this  case,  however,  we  observe  that  the  contribu¬ 
tion  from  the  part  of  the  contour  encircling  the  singularity  at  uj3  is  nonzero.  Therefore  we 
break  up  the  contour  C"  into  two  parts:  the  part  encircling  the  singularity,  G3(r/),  where 
rj  ->  0+  is  the  radius  of  the  circle,  and  the  remaining  part  of  the  contour,  C",  (see  Fig.  2). 

Now,  by  changing  as  before  the  integration  variable,  u  =  icu,  we  can  write 


/  U1  u2  \ 

9&>r)=  eoo0(*>r)  +  |^|/  +  /  j 


du 

2tt 


0  U-+T1 


(U  ~  Uj)  (U  ~  U2) 


u  (u  -  u3 ) 


e  ut  sm 


in  (/(„)  T~) 


—  1 


4nr 


c3(v) 


/d u  ( u  —  Uj^)  ( u  —  u2)  ^ 

27 r  u  (u  —  u3) 

p(u ) 


e-ut  gi  f{u)r/v 


Uj  Ur , 


^+t \h  S  s 


0  u3+rj 


U 


e  ut  sm 


in  (/(«)  r~) 


—  1 


.  r  du 

X«) 

A'k  r  J  2n 

J. 

u 

c3(v) 

e~ut  gi  /(“)  r/v 


(3.28) 


By  changing  the  variable  to  z  =  1  / \/u  —  u3,  we  can  show  that  both  integrals,  over  C3(r]) 
and  C  ,  are  finite  in  the  limit  of  T]  —>■  0+ .  However,  despite  the  fact  that  the  integrals 
are  well  defined,  their  numerical  evaluation  requires  special  care.  We  discuss  details  of  the 
analytical  evaluation  and  the  numerical  implementation  of  Eq.(3.28)  in  the  Supplement  to 
Appendix  Al  of  the  second  yearly  report. 

More  details  on  this  part  of  the  project  were  reported  in  the  second  yearly  progress 
report  and  published  in  [6]. 
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(D.)  Formulations  based  on  the  pulse  temporal  basis  functions  and  band-limited 
temporal  basis  functions 

Temporal  separation  densities  constitute  an  important  element  of  the  full  Galerkin  pro¬ 
cedure  we  use  in  the  evaluation  of  all  of  the  matrix  elements  of  the  components  of  the 
integral  equation  kernel  for  general  dispersive  media.  We  have  constructed  all  pertinent 
temporal  separation  densities  in  the  analytical  form  for  pulse  basis  functions  and  in  terms 
of  one  dimensional  integrals  for  band-limited  basis  functions. 

The  temporal  separation  density  functions 

rT{K) (t)  =  (t1  -  t2)  (3.29) 

depend  only  on  the  relative  time  t  =  t1  —  t2,  and  on  one  relative  temporal  index  m  =  p  —  v, 
where  p  and  v  are  the  indices  of  the  temporal  basis  functions. 

For  the  non-overlapping  pulse  temporal  basis  functions  given  by 


le 


(3.30) 


where  =  p  At,  the  expressions  for  the  temporal  separation  density  functions  (3.30)  are 
quite  simple.  We  find 

-m<1)(‘)  =  ^4  P~*M  -  r-M]  «  3^  *(5) — •  (4r)3’31a) 


r,mV\t)  =  c  J  drrym^(r)  =  H{£)  =  H( 

t 


t  -  tv 
At 


t  -  tr 


ira<4,«  =  I“K)  -  <*(5  - 1)  -  <5(5  +  D] 


[2  8(t  -  tm)  -  S(t  -  tm+1)  -  6(t  -  tTO_i)] 


(cAt)‘ 


dvm{1)(t)  =  =  J_  dvm{1){£) 


dt 


At  d£ 


(3.31b) 

(3.31c) 

(3.31d) 

(3.31e) 

(3.31f) 

(3.31g) 


where 


£  =  ~r - m  =  - (p  —  u) , 

s  A t  At  ’  ’ 

s(£)  =  2Q(£)-Q(£-l)-@(£  +  l), 
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(3.32) 

(3.33) 


(3.34) 


mo  =  (  i-iei)©(i-iei), 


'  0  for  £  >  1  , 

H(  ()  =  ■ 

|  (1  -  £)2  for  0  <  £  <  1  > 

l-|(l  +  £)2  for  —  1  <  £  <  0  , 

(3.35) 

,  1  for  £  <  -1 

^s(£)  =  2<5(£)  -5(£-l)  -8(Z  +  1), 

(3.36) 

§-(m  =  -»(«). 

(3.37) 

{)  =  -MO . 

(3.38) 

In  Fig.  3  we  display  graphs  of  the  the  first  three  temporal  separation  density  functions. 


etas  for  pulse  basis  functions 

Fig.  3:  ?7i ,  772 ,  and  773  functions  for  the  pulse  basis  functions. 


In  the  following  we  discuss  aspects  of  the  derivation  of  the  temporal  separation  density 
functions  (3.29),  in  terms  of  the  band-limited  basis  functions  (APSWFs). 

The  rationale  for  using  such  basis  functions  is  as  follows:  The  block  lower-triangular 
structure  of  the  impedance  matrix  in  time  indices  from  the  causality  of  the  Green  function 
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(G(t,  r)  =  0  for  t  <  0)  provided  the  temporal  basis  functions  (t)  =  T(t  —  tiA t)  do 
not  overlap.  Such  as  restriction  of  the  function  T  to  short  time  support,  however,  is  in 
conflict  with  the  assumption  of  its  band-limited  property,  which  is  essential  in  ensuring 
high-frequency  stability  of  the  solution. 

In  this  situation  the  best  possible  solution  appears  to  make  the  temporal  basis  func¬ 
tions  band-limited  and  approximately  of  finite  (but  not  necessarily  short)  time  support, 
and  restore  the  lower-triangularity  of  the  impedance  matrix  by  other  means  -  namely, 
by  temporal  extrapolation  techniques.  The  required  features  can  be  satisfied  to  a  good 
accuracy  by  the  approximate  prolate  spheroidal  wave  function  (APSWF)  [2]. 

We  assume  that  the  maximum  frequency  in  the  problem  is  /0  =  cn0/(27r).  The  time 
step  is  then  chosen  as 

.17 r 

At  =  o~T  = -  ’ 

2s/o  s  u0 

where  s  ~  10  is  the  oversampling  rate.  This  choice  is  dictated  by  the  condition  that  At 
should  be  of  the  order  of  c  Ar,  where  the  spatial  resolution  Ar  is  of  order  of  Ao/10,  where 
A0  is  the  wavelength  corresponding  to  the  maximum  frequency  /0. 

With  the  parameters  /0  and  s,  the  APSWF  T(t)  is  defined  as 


sin(7rt/At)  sin(o  y/ ( t/N  At)2  -  1) 
7r  t / At  sinh(a)  \J(t/N  At)2  —  1 


(3.40) 


where  N  >  1  is  the  APSWF  width,  and 


a  =  n  N  (3.41) 

is  called  the  time-bandwidth  product  of  the  APSWF.  As  shown  by  Eq.(3.40),  the  APSWF  is 
the  oscillatory  sine  function  modulated  by  a  smooth  factor,  rapidly  vanishing  for  |t|  >  N  At. 

Figs.  4a  and  4b  show  the  APSWF  T(t )  for  /0  =  1  GHz,  s  =  10,  and  N  =  5,  for  which 
At  =  0.05  ns,  and  a  =  14.14.  The  APSWF  is  plotted  there  as  a  function  of  t/ At,  in  the 
linear  and  logarithmic  scales. 
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Fig.  4a:  The  function  T(t)  plotted  with 
t  in  units  of  At. 


Fig.  4b:  The  absolution  value  of  the 
function  T(t)  plotted  in  the  logarithmic 
scale. 


A  difficulty  arising  when  APSWF  are  used  as  temporal  basis  functions  is  that  the  basis 
functions  Tfl  and  Tu  for  different  time  steps  p  ^  v  overlap.  As  a  consequence,  the  resulting 
impedance  matrix  is  not  block  lower-triangular,  i.e.,  A^u  ^  0  for  /j  <  u,  or,  physically, 
that  the  values  of  the  current  at  a  given  times  step  /i  depend  on  current  values  for  future 
times,  v  >  fi.  However,  the  block  lower- triangularity  may  be  recovered  by  extrapolating 
the  current  from  the  past  to  future  values.  A  good  accuracy  can  be  obtained  due  to  the 
fact  that  the  current  is  a  band-limited  function  of  time.  The  result  of  the  extrapolation  is 
that  the  original  (not  block  lower-triangular  matrix)  can  be  replaced  by  a  new  matrix  A1"' 
which  is  block  lower-triangular.  That  matrix  is  computed  as 


(  A^ 


A»u  =  ( 


{  A»v  +  i  AilK  H^-K’K~U 


for  v  <  [i  —  Np  , 

for  fx  —  Np  +  1  <  v  <  /j,  . 


(3.42) 


Here  IPl,u  is  matrix  of  coefficients  obtained  by  solving  a  set  of  least-squares  equations 
which,  physically,  represent  the  condition  that  future  values  of  a  set  of  periodic  band- 
limited  functions  are  represented,  with  the  minimum  error,  as  linear  combinations  of  their 
past  values.  In  Eq.(3.42)  N{  is  the  number  of  required  future  values,  and  Np  the  number 
of  past  values  used  in  the  extrapolation.  While  N{  (typically  about  5)  is  determined  by  the 
width  of  the  temporal  support  of  T(t)  (i.e.,  by  the  APSWF  width),  the  number  Np  of  past 
steps  can  be  adjusted  to  reduce  the  extrapolation  error  to  the  desired  tolerance;  a  typical 
value  of  Np  is  about  10.  Higher  values  of  Np  provide  a  better  accuracy,  but  increase  the 

computational  cost  of  evaluating  the  modified  matrix  (3.42).  On  the  other  hand,  A^u  has 
to  be  computed  only  once,  and  the  resulting  computational  cost  is  not  significant. 


In  Figs.  5a  to  5d  we  display  the  plots  of  the  temporal  separation  densities  computed 
with  the  temporal  basis  functions  (3.40). 
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We  note  that,  since  the  temporal  separation  densities  rj^  have  the  form 


^W)(i)  =  *?°«)(t-A*A0,  3  =  1,2, 3, 4, 

and  the  functions  77°^  (t)  practically  vanish  outside  certain  interval  [— Tmax,  Tmax],  it  is 
sufficient  to  tabulate  them  in  that  interval,  and  construct  other  functions  by  translation. 


Fig.  5a:  The  function  r]°^(t)  plotted  Fig.  5b:  The  function  r?°(2)(i)  plotted 
with  t  in  units  of  At.  with  t  in  units  of  At. 


Fig.  5c:  The  function  r^°(3)  (t)  plotted  Fig.  5d:  The  function  rj°^(t)  plotted 
with  t  in  units  of  At.  with  t  in  units  of  At. 

More  details  on  this  part  of  the  project  were  reported  in  the  third  yearly  progress 
report  and  published  in  [8]. 
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(E.)  Regularization  of  low-frequency  solutions  of  integral  equations. 

We  developed  and  implemented  algorithms  necessary  for  regularization  of  low- 
frequency  solutions  of  integral  equations.  Although  these  methods  are  originally 
formulated  in  frequency  domain,  we  have  been  motivated  by  their  applicability  to  time- 
domain  (TD)  integral  equation  methods. 

One  of  the  critical  problems  in  the  implementation  of  TD  integral  equations  is  their 
stability  for  long  evolution  times.  It  is  known  that  there  are  two  types  of  TD  integral 
equations  instabilities:  those  arising  from  high-  and  low-frequency  phenomena. 

The  high-frequency  instabilities  manifest  themselves  as  non-physical  rapid  oscillations 
in  the  solution,  of  amplitude  increasing  with  time.  They  can  be  reduced  or  eliminated  by 
limiting  the  bandwidth  of  the  incident  signal,  and  ensuring  that  no  high  frequency  signals 
are  introduced  as  the  result  of  time  evolution. 

The  low-frequency  instabilities  are  due  the  appearance  of  low-  (“zero-”)  frequency 
(almost  static)  solenoidal  solutions  which  do  not  generate  tangential  electric  fields,  and  may 
therefore  contaminate  the  correct  solution  (at  least  in  the  formulation  based  on  electric- 
field  equations,  which  are  based  on  the  boundary  conditions  for  tangential  electric  field 
components).  Methods  of  eliminating  these  parasitic  admixtures  are  based  on  separating 
the  solution  space  into  two  subspaces:  (i)  that  spanned  by  divergence- free  (solenoidal)  basis 
functions  and  the  remainder  space,  and  (ii)  the  remainder  space.  An  effective  procedure  is 
then  to  simply  remove  static  solution  components  from  the  subspace  (i)  at  each  evolution 
step. 

Thus,  the  approaches  to  eliminating  low-frequency  instabilities  in  TD  integral  equa¬ 
tions  are  based  on  the  same  solution  space  decomposition  methods  as  those  used  in  the 
low-frequency  limit  in  frequency  domain. 

Our  main  focus  was  on  the  development  of  efficient  numerical  methods  of  analyzing 
topology  of  complex  geometries  and  their  surface  meshes.  These  operations  are  exactly 
the  same  in  the  context  of  frequency-domain  and  time-domain  approaches,  and  in  both 
cases  lead  to  construction  of  of  solenoidal  and  remainder  basis  functions. 

The  main  outcome  of  our  work  was  an  efficient  algorithm  for  constructing  basis  func¬ 
tions  of  the  classes  (i)  and  (ii),  as  mentioned  above.  It  is  applicable  to  triangular  surface 
meshes  and  basis  functions  associated  with  the  mesh  edges.  For  geometries  characterized 
by  topologies  with  a  limited  number  of  connected  components  and  “handles” ,  the  compu¬ 
tational  cost  of  the  construction  is  of  order  0(E),  where  E  is  the  number  of  edges  in  the 
mesh. 

We  have  implemented  the  algorithm  as  a  set  of  geometry  processing  functions  carrying 
out  the  required  topological  analysis  (available  also  as  free-standing  utilities) ,  and  as  a  set 
of  functions  for  constructing  a  “low-frequency  preconditioner”  applicable  to  frequency- 
domain  integral  equations.  The  topological  analysis  routines  can  be  directly  applied  to  the 
analysis  in  time  domain.  Our  implementation  has  been  tested  on  many  complex  and  large 
size  geometries  of  various  topologies,  including  multiple  components  and  many  handles. 

We  believe  that  the  developed  algorithm  methods  are  also  of  interest  in  the  more 
general  context  of  computational  topology,  and  its  applications  in  biology  and  chemistry 
(topology  of  macromolecules,  etc.),  as  well  as  cosmology  (structure  of  galaxies). 
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As  typical  results,  we  show  geometries  created  by  our  topological  analysis  routines 
applied  to  the  surface  (of  genus  4)  displayed  in  Fig.  6.  This  geometry  is  triangulated  with 
about  240,000  edges. 


Fig.  6:  The  geometry  used  in  the  topological  analysis  and  construction  of  basis  functions. 

Fig.  7  exhibits  the  geometry  created  by  our  utility  genloop,  which  constructs  topo¬ 
logically  nontrivial  loop  paths.  These  paths  are  sequences  of  segments  of  the  geometry 
dual  to  the  original  geometry;  they  define  a  set  of  “nontrivial”  solenoidal  basis  functions. 
The  other,  “trivial” ,  loops  are  those  encircling  all  interior  vertices  of  the  mesh. 

Fig.  8  shows  the  geometry  obtained  by  using  the  utility  gentree,  which  creates  tree 
meshes.  These  meshes  specify  basis  functions  belonging  to  the  “remainder”  space  of  non- 
solenoidal  solutions.  In  this  case,  in  order  to  be  able  to  visualize  the  mesh,  we  used  a  much 
coarser  surface  discretization,  with  less  than  4000  edges. 
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Fig.  7:  The  geometry  showing  nontrivial  loops  (defining  solenoidal  basis  functions)  for  the 
geometry  of  Fig.  6. 


Fig.  8:  The  geometry  showing  trees,  defining  “remainder”  basis  functions  for  the  geometry 
of  Fig.  6. 
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(F.)  Fast  Solution  Algorithms  Applicable  to  Finite  Periodic  and  Partially  Peri¬ 
odic  Systems 

We  designed  and  partly  implemented  a  set  of  fast  solution  algorithms  applicable  to  finite 
periodic  and  partially  periodic  systems,  such  as  antenna  arrays. 

The  algorithms  take  advantage  of  the  block- Toeplitz  structure  of  the  impedance  matrix 
resulting  from  the  periodicity  of  the  scatterer.  They  can  be  applied  either  to  the  uncom¬ 
pressed  MoM  impedance  matrix  or  in  conjunction  with  FMM  or  FFT-based  compression, 
resulting  in  what  we  call  Toeplitz-MoM,  Toeplitz-FMM,  and  Toeplitz-FFT  solution  pro¬ 
cedures. 

The  main  effect  of  the  block- Toeplitz  structure  of  the  matrix  is  a  dramatic  reduction 
in  required  memory,  due  to  elimination  of  redundant  matrix  blocks  and  other  storage 
components  associated  with  identical  elements  of  the  periodic  array.  The  matrix-vector 
multiplication  complexity  of  the  algorithms  does  not  change  significantly  (except  that  it 
is  much  reduced  in  the  Toeplitz-MoM  approach;  that  algorithm,  however,  is  in  practice 
applicable  only  to  relatively  small  small  subsystems,  with  <  1000  unknowns). 

We  carried  out  a  thorough  analysis  of  the  storage  requirement  and  matrix-vector 
multiplication  complexity  of  the  three  block-Toeplitz  algorithms,  both  theoretically  (on 
the  basis  of  the  structure  of  the  algorithms)  and  empirically  (based  on  our  preliminary 
algorithm  implementation).  A  typical  result  is  shown  in  Fig.  9:  it  is  the  required  storage 
in  Toeplitz-FMM  and  Toeplitz-FFT  algorithms  as  a  function  of  the  number  of  unknowns 
n  in  a  single  array  element,  for  several  values  of  the  number  m  of  array  elements. 


Fig.  9:  Required  storage  in  Toeplitz-FMM  and  Toeplitz-FFT  algorithms  as  a  function  of 
n  for  several  values  of  m  (corresponding,  e.g.,  to  10  x  10,  30  x  30,  and  50  x  50  arrays). 

These  plots  show  that  the  Toeplitz-FFT  grows  faster  than  for  Toeplitz-FFT  with  the 
increasing  m,  but  slower  (linearly,  rather  than  quadratically)  as  a  function  of  n. 
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Results  of  this  type  and  provide  valuable  information  on  the  regions  of  applicability 
of  the  individual  Toeplitz  methods.  In  particular,  Fig.  10  shows  the  ratio  of  the  required 
storage  for  the  Toeplitz-FMM  algorithm  to  that  of  the  Toeplitz-FFT  algorithm. 


Fig.  10:  Ratio  of  the  required  storage  for  Toeplitz-FMM  and  Toeplitz-FFT  algorithms, 
plotted  as  a  function  of  m  and  n  (in  thousands).  Contours  for  ratios  1  and  5  are  also 
plotted. 

We  can  summarize  the  main  results  by  stating  that: 

-  The  Toeplitz-MoM  approach  is  preferable  for  a  large  number  of  small  elements  (small 
number  of  unknowns  per  element),  where  it  is  much  faster  than  the  other  methods, 
and  the  required  storage  is  still  reasonable. 

-  The  Toeplitz-FMM  algorithm  is  preferable  for  a  large  number  of  “medium-size”  ar¬ 
ray  elements;  a  large  number  of  unknowns  per  element  and  highly  sub- wavelength 
discretization  incur  large  costs  due  to  the  near-field  storage. 

-  The  Toeplitz-FFT  procedure  is  preferable  for  a  moderate  number  of  elements  with 
a  large  number  of  unknowns  per  element  and  possibly  sub-wavelength  discretization, 
approximately  when  the  ratio  n/m  exceeds  10. 

More  details  on  this  part  of  the  projects  were  reported  in  the  second  yearly  progress 
report  and  published  in  [7]. 
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(G.)  Generalization  to  incident  waves  given  as  superpositions  of  Hermite  poly¬ 
nomials 

We  have  enhanced  our  frequency-  and  time-domain  codes  simulation  capabilities  to  in¬ 
clude  numerically  efficient  modeling  of  specific  time-localized  waveforms,  such  as  linear 
combinations  of  Hermite  polynomials  or  approximate  prolate  spheroidal  wave 
functions. 

There  is  a  need  for  such  particular  solver  numerical  capability  in  the  context  of  several 
current  and  future  medical  and  military  potential  applications,  involving,  in  particular, 
detection  of  concealed  small  arms. 

In  order  to  be  in  a  better  position  to  analyze  such  problems,  we  derived  compact 
analytical  expressions  for  the  projections  (on  the  spatial  basis  functions)  of  the  incident 
wave  represented  as  a  linear  superposition  of  Hermite  polynomials  in  temporal  and  spatial 
variables. 

A  particular  example  of  practical  interest  belonging  to  this  class  of  waveforms  is  the 
“Mexican  hat”  waveform,  expressible  as  a  linear  superposition  of  Hermite  polynomials  up 
to  the  second  order. 

We  implemented  the  resulting  formulation  in  the  code  module  generating  the  incident 
wave  projection  on  the  Rao-Glisson-Wilton  (RWG)  basis  functions  defined  on  triangular 
patches,  and  on  pulse  or  band-limited  temporal  basis  functions.  The  RGW  basis  function 
consists  of  a  pair  of  “half  basis  functions”  defined  over  a  triangle  T  and  associated  with 
the  edge  lm  (or,  equivalently,  with  the  vertex  vm  facing  the  edge  lrn) 

'Mr)  =  w-(r  -  rm)  x(r)  ,  (3.43) 

(let 

where  x  is  the  characteristic  function  of  the  triangle. 

The  incident  wave  is  assumed  in  the  form  of  a  linear  superposition  of  Hermite  poly¬ 
nomials  Hrn 

M 

E(inc)(r  ,*)=  JMmHm(7(k0-r-cf))  (3.44) 

m= 0 

of  the  argument  7  (ko  •  r  —  ct),  with  |k0|  =  1,  and  with  a  scale  parameter  7. 

Evaluation  of  the  incident  wave  projection  required  an  efficient  and  numerically  stable 
algorithm  for  computing  convolutions  of  Hermite  polynomials  with  the  basis  functions 
(3.43).  We  have  derived  analytic  expressions  in  the  form  which  is  easily  amenable  to 
numerical  implementation.  Below  we  discuss  the  elements  of  the  derivation  for  the  m  —  0 
term  of  the  series  (3.44).  Each  m-th  term  can  be  constructed  in  terms  of  the  derivative  of 
the  (m  —  l)-th  term. 

The  projection  Bafl  of  the  incident  wave  is  given  by 

Ba „  =  I  as  I  dte-^ko  r-ct)2  M r)M<)  - 

T 
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where  is  a  temporal  basis  function,  chosen  as  a  rectangular  pulse  or  a  band-limited 

prolate  spheroidal  wave  function. 

We  have  carried  out  integrations  over  the  barycentric  coordinates  (£i,£2)  (with  0  < 
Cl  <  1,0  <  £2  <  1,0  <  £1  +  £2  <  1)  analytically.  In  these  coordinates  the  projection  Bail 
is  given  by 

Bafl  =  2A  1 1  dC!dC2  J  dte-^1+{*2+c>2  (a1C1  +  a2£2)$M(<)  , 

T 

where  A  is  the  triangle  area  and  a,  b,  and  c  are  expressible  in  terms  of  the  parameters  of 
the  incident  wave  and  have  quadratic  dependence  on  t. 

We  note  that  the  resulting  expressions  are  difficult  to  program  in  the  double  limit 
when  the  parameters  a  and  b  tend  to  zero  simultaneously. 

We  have  derived  suitable,  numerically  stable  Taylor  expansions  for  these  functions  by 
expressing  them  in  terms  of  auxiliary  functions  Io(a,  b,  c )  and  -Fq (a,  c), 


with 


1  i-4i 

/0(a,  b,  c)  =  J  d£i  J  d£2  erf(a£i  +  6£2  +  c)  , 

0  0 

1  1_?1 

h(a,b,c )  =  J  d£x  J  d£2C1e-«1+^+c)2  =  ^  J0(a,6,c)  , 

0  0 

1  i-«i 

I2(a,b,c)  =  J  d£2  J  d£2e2e-^+^+c)2  =  ^/0(a,6,C)  =  71(6,a,c), 

0  0 


I0{a,b,c) 


F0{a,c ) 


F0(a,c)  -  F0(b,c) 
a  —  b 

H(c  +  a)  —  H(c) 
a 


and 

Jf(Z)  =  y|erf(z)(2z2  +  l)  +  |e-2' . 

For  small  arguments  we  used  the  following  expansions  of  the  functions  H(z)  and  Fo  (a,  c): 


1  1  3  1 


1 


"«  =  i"+r3-60  z>+mzV+ 


F0{a,c) 


^  e  c2  +  cerf(c)  4-  ^  erf(c)  a  +  ~  e  °2  a2 
J2  ce~°2  a3  +  ^QCe~°2  (2c2  -  1)  o4  +  --  -  . 
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Fig.  11:  Representative  display  of  the  current  distribution  on  a  man  with  a  knife. 


Fig.  12:  Comparison  of  the  preliminary  results  for  the  radar  cross  section  for  a  500  MHz 
incident  wave  illuminating  a  man  with  and  without  a  knife 

We  initiated  calculations/analysis  for  representative  scenario  a  man  with  a  concealed 
knife  illuminated  with  different  incident  beams  given  in  terms  of  superposition  of  Hermite 
polynomials  (Fig.  12  shows  a  representative  display  of  the  electric  current  induced  on  such 
a  geometry).  This  analysis  is  in  progress. 
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(H.)  Work  on  aspects  of  foliage  penetration  problems 

We  performed  several  computations  to  verify  the  applicability  of  the  developed  code  to 
the  analysis  of  SAR  returns  from  scenes  containing  objects  of  military  interest  concealed 
under  foliage. 

To  perform  more  realistic  simulations  than  the  plane  wave  (or  superposition  of  plane 
waves)  illumination,  we  added  the  antenna  modeling  option  into  our  integral  equation 
solver.  We  model  the  transmitting  and  receiving  elements  as  delta-gap  structures  built 
into  a  Rao-Wilton-Glisson  basis  function.  The  transmitting  element  at  an  edge,  say,  a 
is  characterized  by  a  voltage  source  Vo  and  a  (connected  in  series)  resistance  Ra.  The 
receiving  delta-gap  at  an  edge,  say,  /3,  is  defined  similarly  to  the  transmitting  one,  but 
without  the  voltage  source,  i.e.,  it  represents  only  the  resistive  load  Rp.  The  goal  of  the 
MoM  antenna  modeling  is  to  obtain  predictions  for  the  antenna  parameters  of  a  single 
antenna  or  a  system  of  antennas,  interacting  with  themselves  and/or  with  other  objects. 
Of  main  interest  are  parameters  associated  with  the  antenna  transmitting  or  receiving 
ports,  such  as  the  input  impedance,  received  voltage,  transfer  function,  etc. 

We  constructed  geometries  composed  of  irregularly  shaped  ground,  trees,  and  a  per¬ 
fectly  conducting  object  located  under  the  trees.  We  also  modeled  the  transmitting  and 
receiving  horn  antennas. 

Below  we  present  some  of  the  relevant  simulations. 

Fig.  13  shows  a  scene  composed  of  a  ground,  four  trees  and  a  perfectly  conduct¬ 
ing  object  of  pyramidal  shape.  The  scene  is  20  m  in  size,  the  trees  are  approxi¬ 
mately  10  m  high,  the  size  of  the  object  is  1  m.  The  computations  were  performed 
at  200  MHz.  The  computational  problem  size  was  approximately  450,000  unknowns. 
The  realistic  material  properties  of  wood,  leaves  and  soil  were  assigned  to  all  the  ob¬ 
jects  in  the  scene.  Figs.  14  and  15  show  the  bistatic  and  the  backscattering  cross 
sections  computed  with  and  without  the  hidden  object.  Fig.  16  shows  the  influence 
of  leaves.  As  anticipated,  this  influence  is  negligible  at  the  frequency  of  200  MHz. 
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Fig.  14:  Bistatic  cross  section  at  200  MHz  corresponding  to  the  scene  of  Fig.  9. 
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Fig.  15:  Backscattering  cross  section  at  200  MHz  corresponding  to  the  scene  of  Fig.  9. 
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Fig.  16:  Effect  of  leaves  in  the  bistatic  cross  section  at  200  MHz  for  the  scene  of  Fig.  13  . 
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Subsequently,  we  constructed  a  scene  composed  of  uneven  ground,  twelve  trees  and 
a  perfectly  conducting  object  of  pyramidal  shape.  The  scene  was  30  m  in  size.  Two 
horn  antennas,  a  transmitter  and  a  receiver  were  placed,  one  on  top  of  each  other,  on 
a  12  m  track  in  front  of  the  scene.  The  computational  size  of  the  problem,  at  200 
MHz,  was  approximately  1,000,000  unknowns.  The  scene  is  shown  in  Figs.  17  and  18. 
Fig.  19  shows  the  backscattering  cross  section  with  and  without  the  concealed  object. 
Figs.  20  and  21  show  the  voltage  drop  at  the  receiver  and  the  transfer  function  for 
the  system  of  the  transmitter /receiver  antennas  computed  at  six  positions  along  the 
12  m  long  track  again,  in  the  presence  and  in  the  absence  of  the  concealed  object. 


side  view  top  view 


Fig.  17:  A  scene  with  a  pyramidal  object  concealed  under  the  foliage  and  two  antennas  (a 
receiver  and  a  transmitter)  mounted  on  a  track. 
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Fig.  18:  The  same  scene  as  in  Fig.  13  with  a  typical  surface  current  distribution. 


backscattering  cross-section,  200  MHz,  hh  polarization 


-30  -20  -10  0  10  20  30 

phi  angle  in  degrees 


Fig.  19:  Backscattering  cross  section  at  200  MHz  corresponding  to  the  scene  of  Fig.  13. 
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Fig.  21:  Transfer  function  computed  at  6  locations  along  the  12  m  long  track. 


The  performed  computations  indicate  the  code  capability  to  compute  routinely  scat¬ 
tering  from  large,  of  the  order  of  one  million  unknowns,  scenes  composed  of  objects  of 
complex  shapes  and  realistic  materials. 

At  present,  the  code  is  being  used  on  an  Air  Force  contract  to  compute  SAR  returns 
and  to  identify  material  properties  of  the  scene  components.  The  scenes  are  approximately 
30  m  in  size.  The  computations  are  being  performed  at  the  frequency  range  from  100  MHz 
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to  1  GHz  and  for  multiple  (approximately  100)  angles  of  incidence.  The  computational 
problem  size  is  approximately  1.5  M  unknowns  at  the  frequency  of  1  GHz.  A  novel  phase 
interpolation  algorithm  (developed  on  the  above  mentioned  contract  and  implemented  into 
the  code)  allows  us  to  perform  computations  every  10  MHz  and  to  accurately  interpolate 
to  1  MHz  intervals.  The  phase  interpolation  algorithm  allows  us  to  synthesize  results  not 
only  for  relatively  narrow-band  pulses  traditionally  used  in  SAR  applications,  but  also  for 
wide-band  pulses. 
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Dr.  Elizabeth  Bleszynski  is  the  principal  investigator  on  this  effort.  Two  senior  researchers, 
Dr.  Marek  Bleszynski  and  Dr.  Thomas  Jaroszewicz  participate  in  the  effort. 


5.  Publications 

1.  E.  Bleszynski,  M.  Bleszynski,  and  T.  Jaroszewicz,  “Fast  Time-Domain  Integral  Equa¬ 
tions  Approach  for  Wide  Band  Pulse  Propagation  in  Dispersive  Media”,  in  Ultra- 
Wideband  Short-Pulse  Electromagnetics,  Kluver  Academic/Plenum  Publishers,  2003. 

2.  E.  Bleszynski,  M.  Bleszynski,  and  T.  Jaroszewicz,  “Fast  Time-Domain  Integral  Equa¬ 
tion  Solver  for  Dispersive  Media”,  in  Ultra-Wideband  Short-Pulse  Electromagnetics , 
Kluver  Academic/Plenum  Publishers,  2006,  in  press. 

3.  E.  Bleszynski,  M.  Bleszynski,  and  T.  Jaroszewicz,  “Block  Toeplitz  Fast  Integral  Equa¬ 
tion  Solver  for  Large  Finite  Periodic  and  Partially  Periodic  Array  Systems”,  IECE 
Trans  Electronics,  Vol.  E87,  No  9,  pp.  1586-1594,  2004. 

4.  E.  Bleszynski,  M.  Bleszynski,  and  T.  Jaroszewicz,  “Fast  Time  Domain  Integral  Equa¬ 
tion  Solver  for  Dispersive  Media  with  Auxiliary  Green  Functions”,  in  proceedings  of 
the  ACES  (Applied  Computational  Electromagnetics  Society)  Conference,  Honolulu, 
HW,  2005. 

5.  E.  Bleszynski,  M.  Bleszynski,  and  T.  Jaroszewicz,  “Application  of  Fast  Time-  And 
Frequency-Domain  Integral  Equation  Methods  To  Simulation  Of  Pulse  Propagation 
Through  Foliage”  in  proceedings  of  the  KIMAS  2003  (IEEE  Conference  on  Knowledge 
Intensive  Multiagent  Systems),  Cambridge,  MA,  October  2003. 

6.  E.  Bleszynski,  M.  Bleszynski,  and  T.  Jaroszewicz,  “Simulation  of  Synthetic  Aperture 
(SAR)  Returns  from  Complex  Scenes”,  in  proceedings  of  the  2005  KIMAS  (IEEE 
Conference  on  Knowledge  Intensive  Multiagent  Systems),  Cambridge,  MA,  May  2005. 


6.  Interactions/Transitions 

Presentations  at  meetings  and  conferences: 

The  following  invited  presentation  were  made  during  the  grant  period: 

1.  “Application  of  Fast  Time-  And  Frequency-Domain  Integral  Equation  Methods  to 
Simulation  of  Pulse  Propagation  Through  Foliage”  IEEE  Conference  on  Knowledge 
Intensive  Multiagent  Systems,  Cambridge,  MA,  October  4,  2003, 

2.  “Fast  Time-Domain  Integral  Equations  Approach  for  Wide  Band  Pulse  Propagation 
in  Dispersive  Media”,  AMEREM  Conference,  Annapolis,  MD  ,  June  2003, 

3.  “Block- Toeplitz  Fast  Integral  Equation  Solver  for  Large  Periodic  And  Non-Periodic 
Finite  Antenna  Array  Systems”,  IEEE  Topical  Conference  on  Wireless  Applications, 
Honolulu,  HW,  October  17,  2003, 

4.  “Application  of  Wavefront-Based  Method  in  High  Frequency  and  Time  Domain  Scat¬ 
tering  and  Propagation  Problems” ,  IEEE  Conference  on  Antennas  and  Propagation, 
Dayton,  OH,  2003, 
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5.  “Simulation  of  Synthetic  Aperture  (SAR)  Returns  from  Complex  Scenes”,  IEEE  Con¬ 
ference  on  Knowledge  Intensive  Multiagent  Systems,  Cambridge,  MA,  May  2005, 

6.  “Fast  Time-Domain  Integral  Equation  Solver  for  Dispersive  Media”,  EUROEM  Con¬ 
ference,  Magdeburg  Germany  July  2005. 

7.  “Fast  Time  Domain  Integral  Equation  Solver  for  Dispersive  Media  with  Auxiliary 
Green  Functions” ,  ACES  (Applied  Computational  Electromagnetics  Society)  Confer¬ 
ence,  Honolulu,  HW,  April  2005, 

Consultative  and  advisory  functions: 

1.  Assistance  to  Air  Force  Research  Labs,  Brooks  Air  Force  base  in  modeling  and  analysis 
of  problems  in  the  area  of  EM  interaction  with  human  tissue  and  foliage  penetration 
(R.  Albanese,  S.  Sherwood,  and  T.  Campbell). 

2.  Assistance  to  Titan  Corporation;  development  of  the  formulation  and  the  code  module 
applicable  to  the  performance  analysis  and  optimization  of  wide  band  antennas  with 
dispersive  elements  (J.  Riordan  and  H.  Sze). 

3.  Assistance  to  University  of  Florida  in  numerical  simulations  of  observables  in  the  con¬ 
text  of  problems  involving  detection  of  tanks  under  trees  (S.  Shabanov  and  T.  Olson). 

4.  An  early  version  of  our  block  Toeplitz  time  domain  integral  equation  solver  has  been 
independently  verified  and  adopted  as  the  Time  Domain  Adaptive  Integral  Equation 
Method  by  the  University  of  Illinois  researches  who  constructed  a  code  based  on  this 
formulation  (A.  Yilmaz,  K.  Aygun,  J.M.  Jin,  and  E.  Michielssen,  “Matching  Criteria 
and  the  Accuracy  of  Time  Domain  Adaptive  Method” ,  in  proceedings  of  the  2002  San 
Antonio  AP  IEEE  Conference,  Vol.  2,  pp.  166-169,  July  2002). 

Transitions: 

The  developed  formulation  and  the  corresponding  solver  provide  a  unique  opportunity 
for  the  simulation  of  such  problems  of  interest  to  the  Air  Force  and  DOD  as  SAR 
scenes  or  identification  of  objects  concealed  under  foliage  or,  e.g.,  on  a  human  body. 
The  developed  technolody  already  resulted  and  was/is  being  used  on  two  Air  Force 
contracts  awarded  to  Monopole  Research. 


7.  New  discoveries 

1.  Formulation  of  a  novel  fast  algorithm  based  on  simultaneous  compression  in  space 
and  time,  and  exhibiting  the  same  competitive  performance  for  dispersive  as  well  as 
non-dispersive  media. 


8.  Honors/ Awards 

No  honors  or  awards  were  received  during  the  grant  period. 
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